#读取tcga数据整理表达矩阵及临床信息文件
#BiocManager::install("BioinformaticsFMRP/TCGAbiolinksGUI.data")
#BiocManager::install("BioinformaticsFMRP/TCGAbiolinks")
setwd("E:\\1.Ovarian\\1.picture\\1.rawdata\\3.TCGA")
#BiocManager::install("TCGAbiolinks")
library(TCGAbiolinks)
library(dplyr)
query <- GDCquery(project = "TCGA-OV",
                  data.category = "Transcriptome Profiling",
                  data.type = "Gene Expression Quantification",
                  workflow.type = "STAR - Counts")
GDCdownload(query)
data <- GDCprepare(query)
save(data,file = "OV_GDC.rdata")

#该data数据中同时包含了多种不同格式的RNA-seq数据，
#可以使用assayNames()函数查看包含的数据名称。
#主要包含以下几种，依次对应index 1-6。
library(SummarizedExperiment)
assayNames(data)
data_tpm <- assay(data, i = 4)
head(data_tpm)[, 1:2]
#data数据对象中还包含了基因注释信息和样本、
#临床信息，可分别通过rowData()、colData()函数提取。
rowdata <- rowData(data) 
coldata <- colData(data)
##提取tpm矩阵
anno=data.frame(id=rowdata@listData$gene_id,
                gene=rowdata@listData$gene_name)
x=merge(anno,data_tpm,by.x=1,by.y=0)
x1=distinct(x,gene,.keep_all = T)
row.names(x1)=x1$gene
tpm.exp=x1[,-c(1:2)]
save(tpm.exp,file = "OV_TPM.rdata")
tpm.exp1=rbind(id=colnames(tpm.exp),tpm.exp)
write.table(tpm.exp1,file = "OV_TPM.txt",sep="\t",quote=F,col.names = F)

##提取临床信息
coldata$days_to_last_follow_up[is.na(coldata$days_to_last_follow_up)] = 0 #is.na()用于返回是否为缺失值
coldata$days_to_death[is.na(coldata$days_to_death)] = 0   
coldata$futime<-ifelse(coldata$vital_status=='Alive',
                       coldata$days_to_last_follow_up,
                       coldata$days_to_death)
c=coldata@listData

phe=data.frame(patient=coldata@listData[["patient"]],
               futime=coldata@listData[["futime"]],
               fustat=coldata@listData[["vital_status"]],
               Age=coldata@listData[["age_at_index"]],
               Gender=coldata@listData[["gender"]],
               Stage=coldata@listData[["figo_stage"]])
phe1=distinct(phe,patient,.keep_all = T)
row.names(phe1)=phe1$patient
phe2=phe1[,-1]
phe2$fustat=ifelse(phe2$fustat=="Alive",0,1)
save(phe2,file = "OV_phebarcode.rdata")
write.csv(phe2,file = "OV_phebarcode.csv")


